%c\delta u-a(x)u+f(u)=0 with zero Dilichlet or Neumann B.C.
%function [J0,g_norm,u,delta_u,solution]=solve
global c a_t b p e t n_t n_p x_t y_t ar nth0 pwr;

nth=input('Input nth Number of the Solution You Want to Find, nth=');
MI=input('Input Minimax Index of the Solution You Want to Find, MI=');
if (MI>1 & nth <= nth0)
    disp('The value of nth is wrong, it will destroy the data saved before');
    nth=input('Re-input the value of nth= ');
    MI=input('Re-input the value of MI=');
end %if

u=zeros(n_t,MI);
delta_u=zeros(n_t,MI);
solution=zeros(n_p,MI);
if MI>1
    [w,delta_w,sol_old]=set_support(nth,MI,bank_u_t,bank_delta_u_t,bank_u_p);
    for i=1:MI-1
        u(:,i)=w(:,i);
        delta_u(:,i)=delta_w(:,i);
        solution(:,i)=sol_old(:,i);
    end
end

[u(:,MI),delta_u(:,MI),solution(:,MI)]=initial;     %Input initial guess

ans=input('Plot the Initial Guess? 1-yes,0-no=');
if ans==1
    %pdesurf(p,t,solution(:,MI));
    pdeplot(p,[],t,'XYData',solution(:,MI),'ZData',solution(:,MI),'ColorBar','off');
end

ans=input('Using Minimax Method? 1-yes,0-no=');

if ans==1
    [J0,g_norm,u(:,MI),delta_u(:,MI),solution(:,MI),u_max]=minimax(MI,u,delta_u,solution);
end

ans=input('Using Newton Method? 1-yes,0-no=');

if ans==1
    [J0,g_norm,u(:,MI),delta_u(:,MI),solution(:,MI),u_max]=newton(u(:,MI),delta_u(:,MI),solution(:,MI));
end

%ans=input('Using Minimax Method? 1-yes,0-no=');
%if ans==1
%[J0,g_norm,u(:,MI),delta_u(:,MI),solution(:,MI),u_max]=minimax(MI,u,delta_u,solution);
%end


bank_u_t(:,nth)=u(:,MI);
bank_delta_u_t(:,nth)=delta_u(:,MI);
bank_u_p(:,nth)=solution(:,MI);

ans=input('Plot the Solution? 1-yes,0-no=');
if ans==1
    resp=input('Generate a spars mesh for graphing? 1-yes,0-no=');
    if resp==1
        [p0,e0,t0]=initmesh(g);
        [p0,e0,t0]=refinemesh(g,p0,e0,t0);
        pdemesh(p0,e0,t0);
        ans=input('Need Another Global Refinement? 1-yes,0-no=');
        if ans==1
            [p0,e0,t0]=refinemesh(g,p0,e0,t0);
        end
        n_p0=size(p0,2);
        for i=1:n_p0
            for j=1:n_p
                err=sqrt((p0(1,i)-p(1,j))^2+(p0(2,i)-p(2,j))^2);
                if err <=0.00001
                    solu(i)=solution(j,MI);
                end
            end
        end
        %     pdemesh(p0,e0,t0,solu);
        %     pdecont(p0,t0,solu,50); hold on; pdemesh(p0,e0,t0,solu-1.5*u_max);
        pdeplot(p0,[],t0,'XYData',solu,'Contour','on','Levels',50,'ColorBar','off');
        hold on;
        pdemesh(p0,e0,t0,solu-1.5*u_max);
    else
        % pdemesh(p,e,t,solu);
        % pdecont(p,t,solution(:,MI),50); hold on; pdemesh(p,e,t,solution(:,MI)-1.5*u_max);
        pdeplot(p,[],t,'XYData',solution(:,MI),'Contour','on','Levels',50,'ColorBar','off');
        hold on;
        pdemesh(p,e,t,solution(:,MI)-1.5*u_max);
    end
end

ans=input('Is the Solution Correct? 1-yes (save), 0-no (no save)= ');
if ans==1
    nth0=nth;
    save  solution_t bank_u_t;
    save  delta_solution_t bank_delta_u_t;
    save  solution_p bank_u_p;
end

